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These lectures notes give an introduction to the fast developing area of research 
dealing with perturbative descriptions of the gravitational instability in an expanding 
universe. I just sketch the outlines of some proofs, and many important contributions 
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(and references) were left out. Many untouched aspects are reviewed in [ p4| which 
also contains a useful reference list. 

In the notes dedicated to the Eulerian approach (§|l| - §[|), I derive the fluid 
equations starting from a microscopic description, I give their linear solution and 
some higher order corrections, and describe a number of applications concerning the 
evolution of the one point probability density function for the density contrast field 
and for the divergence of the corresponding velocity field. 

In the notes dedicated to the Lagrangian point of view 

(§§-©)! recall Zeldovich 

Approximation and its higher orders corrections, and present some applications, like 
their use in describing the statistical effect of the real space-redshift mapping, or as 
approximate description of the full non-linear dynamics. 

1. THE FLUID EQUATIONS 

In this section, I derive the fluid equations as velocity moments of the full microscopic 
evolution equations taken in the mean field, or continuum, limit. We shall restrain to 
the case of perturbations which may be described in the Newtonian approximation. 

In all the following, we use comoving coordinates r = ax, and p = ma 2 9tx, where 
m stands for the particles mass, and a for the metrics scale or expansion factor. In 
that case, the motion and field (Poisson) equations write 

d t p = -V(f>, V 2 = AirGa 2 pS, 

if S = p/~p — 1 stands for the density contrast field. The solution to Poisson equation 
is <f>{x) = —Ga 2 f dx/(p(x') — /0&)/|x' — x|, where p(x') fluctuates in an homogeneous 
and isotropic way around the background value pb- Thus 

-V<j> = Ga 2 J O(x') - Pb ) * ~ * dxt = Go 2 J p(x') * ~ * dxl 

where it is implicit that, in the last integral, one follows the prescription of integrating 
first over angle and only then over x'. 

1.1. Microscopic Description 

Let us denote by /(x, p, t) the one-particle distribution function, i.e. the probability 
density of finding at time t a collisionless particle within the infinitesimal phase space 
volume dp, = dxdp. For N p point particles, each with a trajectory given by Xj(t) and 
Pi (t) , the system is fully described by 

N p 

f(x,p,t) = ^fo(x-Xi(t))tf D (p - Pi(t)), 
1=1 

where 8d stands for Dirac's function. Particle conservation plus Liouville theorem 
(i.e. % = 0) give the equation governing the evolution of the system, % = 0. 
In our discrete particle model, the density is 

p = p b [l + <5(x, £)] = —/ /(x, p, t)p = -3 V 6 D (x - Xj), 

i 

and one recovers the familiar form — V0 = Ga 2 |~z^3 ■ By using the motion 
equation, we then get the Klimontovich equation 

dtf = d t f+-^- V x f - mV0 • V p / = 0, 

ma z 
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a complete description with the field equation 




/. i i > x/ — x 



dx'dp'. 



Macroscopic Description 

Let us now consider a smoothed out (phase space) density obtained by averaging 
over patches of size intermediate between the microscopic level and the system size, 
so that we can take "local" averages. In this "coarse-grain" distribution, we loose 
the discrete character of the particles. This can also be seen as a continuum (or 
thermodynamical) limit such that one formally takes the limit of an infinite number 
of particle, N p — ► oo, while their individual mass goes to zero, m — ► 0, keeping the 
mass density N p m constant. 

In this mean field limit, the Klimontovich equation becomes the Vlasov equation. 
It has formally the same writing, but now the solutions / are to be found in the 
ensemble of smooth, continuous, functions (since the discrete character embodied by 
Dirac function has been lost). Note that by loosing the discrete character of the 
particles, one of course cannot describe the effect of direct encounters (e.g. two- 
relaxation effects). In other words, one only considers the evolution of perturbations 
through fluctuations of the mean field created collectively by many particles. 

Intuitively, this collective description of the evolution should be legitimate when 
the 2-body potential energy at the mean interparticle distance n" 1 / 3 , Grn^n 1 / 3 , is 
much smaller than their average kinetic energy, ma^, which is indeed confirmed by 
more rigorous analysis. 

1.3. Velocity Moments and the Fluid Limit 

A classical method amounts to find equations between various velocity (momentum) 
moments of the one-point distribution function /(x, p, t). These are simply obtained 
by multiplying the equation governing the evolution of / (here the Vlasov equation) 
by products of velocity components(or here momentum components p a ) before 
integrating over all the velocities components. 

The zeroth moment relates the density p(x) = (m/a 3 ) J p/(x, p,t) to the mean 
velocity of a phase space element, v = (p/ma) = J(p/ma)f dp/ J f dp by the 
continuity equation 



which is often rewritten as <9 t <5 + [(1 + 5)v] = 0. 

The first moment obtains by multiplying Vlasov equation by the a component of 
momentum, p a , which yields the velocity equation 



It is then convenient to introduce the centered second moment tensor, i.e. the "pressure 
tensor", P a/3 , defined by 



d t p + 3-p + -V • pv = 0, 



(1) 



a a 



8 t [p a dp + \dp [ P a p p fdp + a 3 pd a 4> = 



(2) 
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to rewrite the source term in the mean velocity equation (g) . 

Truncating the hierarchy: the second moment equation would in turn relate the 
pressure tensor to a "heat flux tensor" (jpp^p 1 ). Of course, one obtains by this 
procedure a hierarchy of equations, each relating the evolution of a moment to the 
one of next higher order. 

Truncating the hierarchy amounts to impose a relation between a given moment 
and the lower ones. In particular, one often assumes an isotropic pressure tensor 

pa/3 p 

= a g — , 

P P 

where 6 a p is Kronecker's symbol, and p is the usual pressure. In that case the first 
moment equation (^) becomes the Euler equation 

^ + V + ( V - V K = _^_M ; 

a a ap a 

which together with the continuity equation and an equation of state relating the 
pressure and the density forms a closed system of equation, often referred to as the 
"fluid equations" . 

Note though that till now we have said nothing about the role of "collisions" . 
And indeed, particles interacting only gravitationally like neutrinos should have an 
isotropic velocity dispersion everywhere, as long as one does not reach the multi- 
stream regime, i.e. before "shell crossing". In fact, early on, the velocity dispersion 
can often be neglected altogether, p — 0, and the equations above describe a so- 
called "dust universe". But of course, one expects on very general grounds (see the 
Lagrangian section for more details) to have an anisotropic, "pancake-like" collapse, 
which cannot be followed any more by these equations (although one may attempts 
to guess an "effective" equation of state mimicking macroscopically the complex 
microscopic dynamics). 

On the other hand, if we have a bona fide fluid made of baryons, their collective 
interactions do maintain a local isotropy and establish an equation of state, with e.g. 
p(p) = p(pb) + c^pbS, where c s stands for the sound speed, c s = (dp/dp) 1 / 2 . In any 
case, the point is that the "fluid equations" can properly be used to describe (some 
of) the evolution of dark matter collisionless particles. 

Other forms of the continuity and Euler equations: it is often convenient to 
introduce the conformal time r defined by dr = dt/a. Denoting the partial derivative 
versus t by an overdot, and redefining v as x, the previous equations become 

j + V[(l + <5)v] =0 and v + -v + (v • V)v = -— - V0. (3) 

a p 

Finally, it is advantageous for some applications to rather introduce D(t) as the time 
coordinate, D(t) standing for the linear growth rate of density fluctuations (see below). 

2. EULERIAN PERTURB ATI VE SOLUTIONS 

2.1. Linear Density Field 

By inserting the first moment equation (|^) into the continuity equation ([!]) one gets 

[pa/3 



d 2 t 5 + 2-dtS = -j V [(1 + 6)V<f>] + -^d a d 
a a z a z 



Pb 



+ {1 + S)v a v 13 



(4) 
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The Eulerian perturbative regime is when 

/ .\ 2 

and (— J ~5<1, 

where d stands for the coherence length of the perturbation and t its dynamical time, 
t ~ (Gp) -1 / 2 (in an fl = 1 universe, the FRW equations lead to 6nGpt 2 = 1 in the 
matter era). The second requirement is simply that gradients should be small. In 
linear perturbation theory, one neglects all term in 5 2 ,v 2 ,5v 1 and the fluid equations 
assume the form 

dtS + V • va = and d 2 S + 2-d t 5 = ^4 + ^Gp b S. (5) 

a p b a 2 

Solutions are known for many different cases, i.e. for different equation of state, and 
different time evolution of the factor of expansion a. Here we recall some of the most 
useful ones. 



2.1.1. Dust Universe, p = 0: since there are now no spatial derivatives in 
the equations, the evolution is self-similar for arbitrary density contrasts 5. Let us 
denote by D a and D b the two growth rates solution of the second order differential 
equation, and by A and B the associated spatial parts, so that the solution if the 
linear combination 

<f(x,t) = A(x)D a (t) + B(x)D b (t), 

D a corresponding by convention to the fastest growing mode 

In the Einstcin-dc Sitter case, f2 = 1 and a oc t 2 / 3 . Thus eq. (Q) becomes 

d 2 8 + -^dtS = jjfj and we have D a — t 2 ^ 3 and D b = t~ x . Note that the exponential 
growth of a static universe (when d = 0) is now much slower. 

In the open case, full solution may be obtained by changing the time coordinate 
to a given by 



8irGpb 



1 + (flo 1 - 1) - 



Of) 



But in the asymptotically late regime when the curvature term dominates, i.e. a act, 

}d t 8 = 0, since 4irGp b < r 2 . Thus D a = 1 



one simply has d 2 5 + \d t S + 



d 2 t 5- 



and D b = t ] At best the density contrast is frozen in the free expansion regime. 
And this happens when (fi^ 1 — l)a/a,Q ~ 1, i.e when 1 + zj reeze — fi^ 1 — 1. In other 
words, a perturbation grows like t 2 ' 3 at z > Zf reeze and stops afterwards. 

Fixing j4(x) and -B(x): the solution to our second order differential equation is 
fully specified by giving, e.g. the initial density contrast Si and velocity Vi at time U, 



AD a (ti) + BD b (ti), and 



V'V, = . dD a 



a, 



dt 



B 



dDb 
dt 



Introducing the Wronskian W = D a {U) ^§r\ - - D b (t t ) ^§r\-, we have then 



5 = 



W 



d_Dt 

dt 



D h 



dDa 

dt 



[D a D h {U) - D b D a (ti)} 



() 



In the Einstein-de Sitter case, W — — | i?g-W2iM*i) , if the initial velocity may be 
neglected, one then finds 

\D a + \D b 

with 8 = ^D a Si at late times (in the following, we shall drop the subscript a, and D 
will refer to the fastest growth rate). On the other hand, in numerical simulations, 
one usually sets the velocity field to select the growing mode only. 

2.1.2. Jeans length of an ideal fluid withp = p(p): plugging p(p) — p(pi J ) J rc 2 pi ) 5 
in (JsJ) , we get 
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d 2 5 + 2-dtS = f -) V 2 S + AnGp b 5. 
a V a J 



a \ a . 

It's Fourier transform is 

2 



a 



d 2 t 5 k + 2-d t S k 



4nGp b 



fa, 



which shows that there is a characteristic wavelength Aj = 2n/kj = ir-j^ ~ c s t for the 
scale of a perturbation. Above this Jeans length, the pressure cannot counteract the 
self-gravity of the perturbation (c s t < A), while at A < Aj the perturbation oscillates 
like a sound wave. 

2.2. Linear Peculiar Velocity Field 

2.2.1. Solution in a dust Universe: the linearized continuity and velocity 
equation in the pressureless case reduce to 

d t 5 + —Vv — and <9 f v + -v = — Y^t_ 
a a a 

Since S oc D, we have V • v = ad t 8 = —HafS = — 47r g°^ a 2 V 2 </> where H stands for 
Hubble's "constant", a/a, and /(f2) = (a/D) dD/da can be well fitted by fl 6 over a 
large range. Thus 

v _ 2 / -V0 V x J7(x) 
3 -ff 51 a a 

which satisfies the first equation for an arbitrary t/(x). This last term decays as a -1 , 
as does v in the absence of a driving term. Since it quickly becomes negligible, we 
find that at the linear order the velocity is parallel to the acceleration at late times, 
the proportionality constant depending rather strongly on the density parameter. 

2.2.2. Application to data: the strong Q dependence of f(£l) has been used in 
the context of the linear bias model, in which one assumes that, at least on sufficiently 
large scales, the density contrast field computed from galaxies, S gi is simply on average 
proportional to the underlying density field, i.e. S g — bS (although that may turn out 
to be a poor approximation for some applications). In that case, the velocity field is 
proportional to the parameter = f(Q)/b. Since one may estimate 8 g through 

i r 1 V- 1 c / x . . i 1 V- 1 

1 + d g = — ^7 — -°d\x- — Xj) with n = — 
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where the density n = y ^/^.a obtains from the catalog selection function <fr, one 
may then compare the predicted velocity from the observed galaxy density field, v pre d 
to the "observed" one. One thus try to find the f3 which best match the two velocity 
fields, with 



_ f3Ha 

Vpred - 



f 5g(x ')W/(|x'-x|)i^riV. 



Note that in practice, we do not have an infinite catalog, and we have to truncate 
the integral at a maximal radius x ma x selected such that the signal still dominates 
over the shot noise (i.e. beware of 1/$ !). Even worse, not all the sky is surveyed 
homogeneously and isotropically (e.g. close to the galactic plane), so that the density 
field has to be extrapolated in some way. Note also the introduction of a window 
function W (whose size is typically greater than 500 km/s) to average over a rather 
large volume in order to reduce the level of non-linearities, reduce the shot noise, the 
distance uncertainties and the impact of the triple-value zone in relating redshifts to 
distances (ax — cz — v is solved iteratively) 

Alternatively, one could use the proper velocities of galaxies to compare with the 
observed density field. This area has been the subject of intense research in the last 
few year, and will be covered in much greater detail in the notes of A. Dckcl (this 
volume). For a more in-depth view, one could look at the proceedings |l(| of the 
"Cosmic Velocity Field" meeting held in 1993, or the recent review [p7| . 

2.3. Non-Linear Corrections 

Let us now come back to the density field equation (^) for a pressureless medium, 

d 2 5 + 2-d t 5 = -iy [(1 + tf)V0] + \d a dp [(1 + 6)v a vP] . (6) 

In order to find a description of the evolution which may remain valid longer than 
the linear one, we write the density contrast and the velocity fields as perturbative 
expansions 

<f(x,t) = <5 (1) (x,t)+<5^(x,t) + ... , 

v(x,<)=v«(x,i)+v( 2 >(x,i) + ... , 

where S^(x,t) and v^'(x,t) are the linear terms obtained previously. They are of 
order S t = e (if we write $W(x,t) = s(x)D(t), vW a (x,t)) = (f/H)d a <$> h if stands 
for the initial potential J > an d we retain only the fastest growing modes). 

The next term 8^ = 0(5^) 2 — 0(e 2 ) is a non-linear correction, solution of the 
equation when only quadratic non-linearities are retained. The source terms are then 
products of the linear solutions. The resulting second order solution can then be used 
iteratively to obtain still higher order solutions... 
The solution to the equation of order two in e is 



<S^ = D{tf 



^[1 + K(t)]e 2 (x)+Ve(x)-VA 



T 



(7) 



where T 2 = J2 a 0=i{T~uf3) 2 is the contraction of the (linear order) tide tensor 



r aP = VoV^, - ^5 a0 V 2 <^ l (8) 
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whose presence can only speed up the collapse. The parameter n(t) is a slowly varying 
function of cosmological time which is well approximated (Bouchet et al 1992) by 



14 



(9) 



in the range 0.05 < < 3 (the accuracy of this approximation is then better than 
0.4%). For Q = 0, k — j. The exact expression for n(Q), valid in the entire range 
fl > 0, is given in ||). For fl = 1, k = and we find 



the classical Einstein-de Sitter solution (e. 



hereafter LSS, eq. [18.8]). 



Vorticity and Shear By taking the curl of Euler equation (|3|) when the conformal 
time r is used, we find that gravity alone cannot generate vorticity. Indeed, introducing 
lu = V x v, we get (an overdot denoting partial derivatives versus r) 

_ , x a Vp x Vp 

LU = V X (V X LU) LU . 

a p A 

Thus lu = for a pressureless fluid which is initially vorticity-free. Furthermore, the 
density field equation (ra) for a pressureless medium can be cast in the following form 



4 S 



3 1 + S 



Q + 



a a pa 



af3 



1 



-LU 



inGpbS 



T af3 _ 1 
2 



0x Q 



+ 



Thus the 



where the velocity shear tensor now has the form a € 

shear accelerates the collapse, while vorticity slows it down. This is this acceleration, 
absent at the linear order, that showed up at the second order [eq. (0)] through the 
(linear) tide tensor [T af) = -^ct q/3 - % i ^ x )- Finally, let us stress that the rate of 
collapse depends on the three eigenvalues of the tide tensor, which implies that the 
highest peaks of the density field do not compulsorily collapse first (see §5.3). 



3. SOME BASIC APPLICATIONS 



3.1. The PDF as a Measure of Large Scale Structures 

Many approaches can be used to characterize statistically the observed large scale 
structures in the Universe. Here we focus on the one-point distribution function, 
or PDF. Empirically, this approach was first applied by Hubble in 1934 Let 
Q(r) be a physical quantity whose value depends on the location. Let Qi(r) be the 
smoothed value of Q on scale £, i.e., Qe{r) is obtained by convolving the scalar field 
Q(r) by a smoothing window Wt of characteristic scale I {Qi{r) — Q(s) * We(r — s) = 
/ d 3 sQ(s)Wi(r — s)). The probability that Qe lies between Q and Q + dQ is 
P(Q, I) dQ, P being its PDF. 

If Q = N where N is a number of galaxies, and W is a top-hat window, then the 
PDF gives the probability of finding N galaxies in a sphere (in 3D) of size I. In the 
following, Q will be either the mass density contrast field S = p/~p—l, or the divergence 
of the associated velocity field in units of the Hubble constant, 9 = V • v/_ff; we will 
consider the two most commonly used types of windows, a top-hat and a gaussian one. 
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Later on, we shall discuss the relation of the mass distribution to the distribution of 
luminous galaxies. 

It is frequently assumed that, very early in the history of the Universe, the density 
contrast s = 6(t = i,-) could be taken as normally distributed, i.e., 

-2 

P(e, £) oc exp 



where of (£) = (ef) stands for the initial variance of the density contrast field smoothed 
on scale I. It is assumed to be small at all scales, ct, <C 1. The second moment is the 
only quantity we need to know in order to fully characterize such a field. Indeed all 
odd moments of a gaussian PDF are zero, while the 2n-th moments are proportional 
to the n-th power of the variance, (e 2 ") oc of 1 (the corresponding specification in 



Fourier space is given in (13)). In the following, we shall refer to such a state as 
gaussian initial conditions. 

Our goal now will be to derive the properties of the PDF today, under the 
assumption of a time-evolution under the sole influence of the gravitational instability 
acting on gaussian initial conditions. Under the influence of gravity, underdense 
regions become even more underdense (although not indefinitely, since the density 
is bounded from below), while positive density enhancements tend to grow without 
bound. Clearly the symmetry of the distribution cannot be maintained, and the PDF 
becomes skewed, i.e., (<5 3 ) departs from zero. The distribution also develops a non-zero 

kurtosis (i5 4 ) — 3 (S 2 ) 2 . In the limit of a small variance, (5 2 ) <C 1, the development of 
skcwness and kurtosis are the most important effects. We shall concentrate on those 
first. 

3.2. Moments of the PDF 

3.2.1. Skewness of the Unsmoothed Density field: by using the perturbative 
results in the weakly-nonlinear regime, we can calculate the gravitationally induced 
skewness under the assumption that e is a random gaussian field. The lowest order 
terms in the series for (<5 3 ) are 

(S 3 ) = ( x S^ 3 ) + ^d^ 2 S^)+0(s 5 ) . (10) 

The linear solution implies that the first term is D(t) times the initial skewness, 
which is zero for gaussian initial conditions. Similarly, at linear order, all higher order 
moments will be proportional the initial ones; thus a normal PDF remains normal in 
linear theory. The second order solution (0) shows that the second term involves (e 4 ), 
which is cx cr 4 for gaussian initial conditions. Thus the skewness ratio 

S 3 = ±—L ( n ) 

is a constant versus scale, and one finds (ref. ||) up to terms of order e 2 (since 
<£ 5 ) = 0), 

S3 = 4 + 4^)^y + ^ 2/63 -l) ■ (12) 

The ~ sign above applies to the range of applicability of equation (||) . The first term of 
this equation, 34/7, had been obtained by Peebles more than a decade ago (LSS, §18). 
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The weak f2-dependence of the full expression shows that nearly all the ^-dependence 
of the skewness (S 3 ') comes from that of the square of the variance. It simply reflects 
the fact that the second order growth rate is nearly equal to D(t) 2 , as can be seen 
from (0). 

3.2.2. Smoothing the density field: to make contact with observables, we want 
the skewness of the density field 5i, when smoothed with either a top-hat or a gaussian 
(spherically symmetric) window, which satisfies 

. 2 j3„ _ if! 



W(x) (Fx = 1, and J W(x)x 2 d A x = i 

The first equation insures a proper normalization to unity, while the second requires 
the effective half-width I to be finite. The top hat case is appropriate for comparisons 
with the observed frequency distribution of galaxies. (In that case, discreteness 
corrections must be taken into account before comparing with the theory, e.g., 
((N/N- l) 2 ) = 1/N+ (S 2 ) if one adopts the Poisson model, see e.g., LSS §33. 
A gaussian window, on the other hand, removes these discreteness fluctuations, since 
it does not have sharp edges.) Here we recall the main results of |3lj ]. 

The calculations are most conveniently done in Fourier space. For the initial field 
e, the Fourier components are given by 



/d 3 k d 3 k' 
—— A P(k) P(k') W k W k , WV k1 T(k, k') + 0(a 2 ) . (14) 



For gaussian initial conditions, one has 

(£ k e k '} = S D (k + k')P(k), <£ k £ k '£ q ) = 0, 

(13) 

<£k£k'£ q £q'} = P(k)P(q)6 D (q + q')6 D (k + k') + cycl. (2 terms) , 

where 6n is the Dirac delta, and P(k) is the (initial) power spectrum. Now we can 
use jib] ) and (|t|) to derive S3 to lowest no n- vanishing order, 

d 3 kd 3 k' 

Wf 

Here T(k,k') stands for T(k,k') = 4 + 4«(fi) - 6ju(fc/fc') + [2 - 4«(f2)] P 2 ((J>), where 
/1 = k • k'/fcfc', and P2 is a Legendre polynomial. If there is no smoothing (Wk = 1), 
the dipole and quadrupole terms integrate to zero, and one simply recovers (|l2|). On 
the other hand, as soon as one introduces smoothing, the result does depend on the 
initial power spectrum P(k). 

Let us assume that P(k) oc k n . Then, for a top-hat smoothing, and —3 < n < 1, 
the equation ( |l4"| ) yields after painful calculations the simple result 

S 3 = 4 + 4/e(fi) - (3 + n) . (15) 

Actually, a careful inspection of the expression ( fl4| ) shows that S3 should only 
depends on the effective (logarithmic) slope of the power spectrum at the smoothing 
scale. This was confirmed in |3l| (see also figure [l]) by comparing numerical integration 
for a CDM power spectrum with a prediction using @. The result above was 
generalized in p| to an arbitrary gaussian field, smoothed with a top hat filter on 
a scale £, 

S 3 (£)=4 + 4 K (n)- le , with 7£ = ° x / . 
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For a pure power law P(k), we have 7 = 3 + n, in agreement with (|l5|). 

The rather messy calculations involved in computing the case of smoothing by a 
gaussian were recently completed in J3(| who express it in terms of hypergeometric 
functions: 

S3 = 3 2 Fi \ m, m, ^, ^ - + 2-P1 fm, m, |, (16) 
with m = 2i 5^. The behavior of this expression with the spectral indices n and £1 is 



in fact very similar to that of the top-hat case, see figure 7.2 below 



3.2.3. Skewness of the Velocity Field PDF The continuity equation gives the 
divergence of the velocity field in terms of the time derivative of the density contrast. 
Let us call T3 the skewness ratio for the divergence of the velocity field in units of the 
Hubble constant (9 = V • v/H). One then finds || for a top-hat smoothing, 



r„ = if) = _ -°- 6 



5)2)2 



26 

— - le 



(17) 



The corresponding (complicated) expression for the a gaussian window has recently 
been obtained (together with the corresponding S3) in [j36f . 

The asymmetry in the distribution of 9 is directly related to the asymmetry in 5: 
voids and clusters in the mass distribution correspond to sources (9 > 0) and sinks 
(9 < 0) in the velocity field. This quantity has the same dependence in fl than the 
ratio of the predicted versus observed velocity in density-velocity field comparisons, 
but it uses an intrinsic property of the field, thereby bypassing reference to the density 
field and to an unknown biasing parameter. 

It turns out that 9 may be advantageously measured by using Voronoi or Delaunay 
tesselations of space |t| . It remains to be seen whether these methods provide accurate 
and unbiased estimators of T3 when the number of sampling points is comparable to 
those of observed catalogs, and one has only the line-of-sight velocity (i.e., bypassing 
the use of a POTENT- like machinery). First analyses |2i| are encouraging in that 
regard. Finally, let us note that an analytical expression for T3 in the gaussian 
smoothing case has recently been obtained (together with the corresponding S3) in 



3.2.4- Kurtosis Using similar techniques to those described above, one can 
compute the Kurtosis coefficient 



(S 4 ) - 3a 4 



and the corresponding T4 for the expansion scalar 9. One finds 
smoothing 
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This was actually derived an equation for the evolution of the full generating function 
(see below). The corresponding expressions for a Gaussian smoothing can be found 

in m. 
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3.3. Edgeworth Expansion of the PDF 

As was shown in [^2| , the previous results on the moments of the PDF in the weakly 
non-linear regime can be used to examine how gravitational instability drives a PDF 
away from its initial state. Let us start with the density field PDF, P(5,£), or rather 
p{v), the PDF of the density field in terms of the standardized random variable 
v = Se/cre. Since we want to describe an evolution from gaussian initial conditions, it 
makes sense to consider an expansion of p(y) in terms of <j){y) = (2n)^ 1 / 2 exp(— z/ 2 /2) 
and its derivatives. The Gram-Charlier series (Cramer 1946 and references therein) 
provides such an expansion: 

P {v) = co 4>{v) + ^ 4> (1) iy) + | ^ (»/) + ..., (18) 

where c m are constant coefficients. Superscripts denote derivatives with respect to v. 
d m <h 

4h\v) ee |-| = (-irH m (u)^(u), (19) 

where H m is the Hermite polynomial of degree to. The Hermite polynomials satisfy or- 
thogonality relations (e.g. Abramowitz & Stegun 1964), J_ H m (v) H p (u) <j>{v) dv = 
0, if m 7^ p, and = m ! otherwise. Therefore, multiplying both sides of equation ([l8]) 
by H rn and integrating term by term yields 

p oo 

Cm={-l) m / H m {v)p{v)dv. (20) 



Equation (|20j) gives Co = 1, c\ = C2 = 0, while for the next four coefficients in the 
series we obtain 

^ = (-1)™ S m crf - 2 , for3<m<5; cq = S 6 aj + lOSf erf . (21) 

Thus the S m {£) have both a dynamic and a static application: they describe the time 
evolution of moments of the PDF at a fixed smoothing scale £, and they also describe 
the relation between moments of the PDF at a fixed time on different smoothing 
scales. In the latter case, one must also include the scale-dependence of the S m if the 
initial power spectrum is not scale-free. 

We have seen that perturbation theory and gaussian initial conditions imply that 
(<5|) oc er^, and S3 is therefore an "order unity" quantity when o~e is the "small" 
parameter. The same is true for all remaining reduced moments, S m — 0(1) for all 
to j| & ]22| . It implies that the Gram-Charlier series is not a proper asymptotic 
expansion for p{v). In an asymptotic expansion, the remainder term should be of 
higher order than the last term retained. However, if we truncated the series (|I|) at 
the term, which is 0(c 2 ), we would miss another 0(o~ 2 ) contribution coming from 
cq (equation (|2l|)). In order to deal with this problem, let us rearrange the Gram- 
Charlier expansion by collecting all terms with the same powers of a. The result is 
the so-called Edgeworth series, with the first terms given by 



pM = ^)-|j^ (3) W- 



-0(a 3 ) .(22) 



Cramer |L7|) lists the Edgeworth series to higher order, and he proves that it is a 
proper asymptotic expansion. This proof is directly relevant to our purposes, since it 
implies that there are no additional 0(a 2 ) terms hiding in the Gram-Charlier series at 
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m > 6. The Edgeworth series thus provides a series expansion for the evolving PDF 
in powers of the r.m.s. fluctuation a. 

Of course, if the variable we are interested in is 8 instead of S, one just adopts 
the appropriate analogues of S m , e.g., S 3 — ► T 3 = (# 3 ) (# 2 ) , and so on. Similarly, 
if we are interested in the PDF shape in rcdshift space, one just needs to use the 



corresponding analogue of S m , S^, see §7.2). The Edgeworth series may also be used 



to relate S3 to other measures of asymmetry like (S\S\) which, according to J4l| may 
offer better signal to noise ratio than <^(5 3 ) when applied to real galaxy surveys (by 
being less sensitive to the tail of the PDF, i.e., to rare events). This latter quantity is 
not easy to compute directly by perturbation theory (see the appendix of (3^] ) , but it 
is trivial to obtain 

(6\5\) = \hrS3<T 3 + 0(a 5 ). (23) 
v »7r 

by using the Edgeworth series. 



3.4- Validity of Perturbation Theory 

The previous results assume that the system never gets to be strongly non-linear, i.e., 
an <C 1 is true at every scale. But in the real universe, we know that the observed 
density contrast is very large at small scales. One only gets (T(<<1 for I 8/i _1 Mpc. 
It is thus by no means obvious that there is any range today for which perturbation 
theory might be applicable. On the other hand, it has long ago been noticed that 
linear perturbation theory yields a good description of the variance [or the 2-body 
correlation function oc D(t) 2 / D{ti) £(U)] measured in fully non-linear numerical 
simulation. 
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Figure 1. The measured Ss(£) in a CDM simulation (triangles) is compared to the 
theory (circles; the error bars come form the numerical uncertainty in evaluating the 
integrals giving S3). The left panel corresponds to a spherical top-hat smoothing, 
while the right one corresponds to a gaussian smoothing. The vertical dashes mark 
the limit between the strongly non-linear regime (a > 1) and the weakly non-linear 
one (cr < 1). Courtesy Colombi 1993 & Juszkiewicz ei al 1993a. 



It turns out that what holds true for linear perturbation theory results is also true 
at higher orders. Indeed, the top-hat expressions for S3 have been checked against 
the simulation results in [|l9| , Q & Q , in the case of scale- free initial spectra. In the 
case of a gaussian smoothing, the theory was checked with the scale free-simulations 
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in pgfl . In the CDM case, comparison were made with the results in as well as 
simulations performed for that purpose, see figure [l]. In all those cases, the agreement 
between the perturbation theory and N-body experiments was excellent (sec |3l| for 
further details), up to surprisingly large values of the variance ~ 1. 









,1111,1111,1- 


-1 




^^^^ ~: 


-2 


- 1 1 : 






:/ | : 




-3 


f , ; 




-4 






p. 







o 



is 



o 



M 
O 



10 




.1,1111,1111,1 


1 1 - 




















1 ' \: 




: 1 i 











10 



Figure 2. Comparisons of the Edgeworth PDF with simulations (solid linos), 
details in the text. Courtesy Juszkicwicz et al 1993b. 



Sec 



One can also check when the Edgeworth series and perturbation theory are 
appropriate to describe the overall shape of the PDF. The solid lines of figure || show 
the PDF of the smoothed density fields of matter density (left column) and velocity 
divergence (right column), at three different stages in scale- free n=-l simulations. The 
corresponding rms amplitude (<r and ag) are shown in each panel, and the variables 
on the horizontal axes are v = 8/a and v = 9/ag. In all the panels, dotted and dashed 
lines show the approximations keeping the terms of order a and a 1 respectively in 
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eq. (|2^)- The approximations work well for S^a < 1 and \vo~\ < 1 (and equivalent 
requirements for 9). They begin to break down outside of that range, as expected. 

So far, all the comparisons done lead to the same conclusion: perturbation theory 
predicts remarquably accurately the simulation results when a £ 1 i\va\ < 1 for the 
Edgeworth PDF), even when small scales are strongly non- linear. All the discrepancies 
found could actually be traced to limitations of the simulations themselves. What 
remains to be found though, is why perturbation works so well, even at scales where 
the density contrast is comparable to unity. 



4. FURTHER DEVELOPMENTS 

4-.1. Higher Order Moments 

While appropriate to break fresh grounds, the pedestrian calculation approach used 
so far would be quite unwieldy to compute higher order moments. Fortunately, 
Bernardeau || found an elegant shortcut by means of a dynamical equation for 
the generating function of quantities closely related to the moments. This function 
simply obeys the evolution equation of the spherical collapse. One nice aspect of this 
approach is that it also gives the moments hierarchy for approximate descriptions of 



the dynamics (like those described in § 7.3 and afterwards). Indeed, in that case 



the moments are then obtained from the generating function obeying the evolution 
equation of the spherical collapse, as given by the considered approximation ]io| . 
In addition, this approach was successfully generalized to account for a top-hat 
smoothing [Q. Unfortunately, this generalization uses a trick specific to the top- hat 
case. Gaussian smoothing calculation still have to be done "by hand" . used for 

4-- 2. Higher Order Corrections to the Moments 

All the calculations presented so far were performed retaining only the leading non- 
linear corrections, that is in the limit a — > 0. Recently, though, the next higher order 
correction to all moments was obtained jl6| , albeit in the absence of smoothing. These 
(one-loop) corrections start to dominate the leading non-linear (tree-level) corrections 
when a ~ 1/2 for the density field and a ~ 1 for the divergence of the velocity field. 

The only case so far when these one-loop correction were obtained in the presence of 
smoothing is for the simplest moment, the variance itself ]37| . One can then approach 
the long-standing "previrialisation" problem: do the presence of small scale slows 
down the collapse of larger structures encompassing them? Or in other words, do 
non-linear interactions between different scales increase or decrease the rate of growth 
of structures. Perturbation theory can at least give us a hint of what the first non- 
linear corrections do. Lokas et al |37| find that the weakly non linear growth of 
variance is very similar to the linear one, when the index n of the primordial power 
spectrum is close to -1, in agreement with the numerical simulation results of [^| 
who concluded to the absence of a previrialisation effect. But it is decreased if it is 
shallower, i.e. if small scale structures are more prominent. This does bring support to 
Peebles' (very educated) guess when he argued that small scale structures will promote 
the development of non-radial motions in the larger collapsing region, these motions 
acting effectively as a stabilizing pressure term. Once again, those perturbative results 
are confirmed by simulations results, as can be seen on figure |[ 
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4-3. General Initial Conditions 

50 far, we focused our attention to the case of Gaussian initial conditions for which 
the reduced moments, or cumulants, are initially zero, and thus S3 and S4 too. Wc 
found that in this specific case, the 5*3 and S4 which develop under the influence of 
gravity take on particular constant values. But what happens for more generic initial 
conditions? Can the predictions for the Gaussian case can be used as observational 
tests of the Gaussian hypothesis? 

Initially, one expects on dimensional grounds to have ((5 3 ) oc a 3 , i.e. S3 oc l/o~. But 
non-gaussian fluctuations will also built a gravitationally induced contribution during 
their evolution, which may quickly dominate the initial term. This was demonstrated 
in specific cases by |l5| & and |24| succeeded in performing the calculation for 
arbitrary initial conditions. They found that s = (5 3 ) /er 3 is equal to its initial value, 

51 plus a term proportional to er, the coefficient of proportionality depending on 
the assumed initial conditions. Unfortunately, biasing (see below) may change this 
proportionality constant, thus comparisons to data will only put an upper limit on Sj. 

Fortunately, as shown in similar calculation at the next order show that 

k = (<5 4 ) /cr 4 is equal to its initial value ki plus a term linear in a. This term is zero 
in the gaussian case (the first non-zero term is in a 2 , leading to a constant £4). This 
measure should thus be a more powerful probe of the gaussian hypothesis, since then 
this is the scaling with a which is altered in the generic case, not a coefficient in front. 
The gravitational evolution is then less likely to cast a shadow on the signature of 
non-gaussian initial conditions. 



4--4- Biasing and comparisons to data 

As was mentioned earlier, it might be that galaxies are not fair tracers of the mass, 
even at the large scales where perturbation theory is applicable. But let us suppose 
that this biasing is local, i.e., that the galaxy field at scale £, S g j is some unknown 
function B of the underlying field at the same location, 5 g ^(x) = B[S((x)]. Then a 
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Taylor expansion would give 



* 9 = E 

k=0 



-6 k 

k\ ' 



where bo is set so that (5 g ) = 0, b\ is the linear bias introduced previously, and 
bk = d n B/dy n \ y= o. Replacing this into the expression for skewness, one sees that the 
galactic field skewness, 5*1 , is given by [B2| 



Si 



S3 
h 



3h 



(24) 



If the biasing is not scale-dependent (why should it be on very large scales?), then the 
constancy with scale of S3 is preserved. This agrees well with the determination flll"| 
in the IRAS 1.2Jy catalog (see fig. |]) or in the SSRS |p6[ . In addition, if a linear 
bias is assumed, like it is customary in density-velocity comparisons for measuring 
(3 = r2 6 /6i , then one gets the linear bias parameter fei, independently of the (poorly 
known!) value of Q. Indeed, once the shape of the window function is fixed, S3 
depends in a known way on the local slope of the power spectrum only, which can be 
determined self-consistently by measuring the slope of the second moment in the same 
catalog. 
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Figure 4. Skewness measurement in the IRAS 1.2 Jy catalog. Courtesy Bouchet et 
al 1993) 



Of course, there might be a non- linear biasing term, 62, and a skewness 
measurement alone will then only constrain the combination (^4j). One cannot 
hope to break this degeneracy by measuring the galactic field kurtosis, since then 
a new unknown, 63, would enter the game. On the other hand, one can compare 
measurements in different catalogs, e.g. optically and infrared selected, to obtain ratios 
of biasing parameters |23 (see also pal). In addition, by extending systematically 
considerations similar to those above, p3fl develop a general (local) bias transformation 
theory and show that the signature of gaussian initial conditions evolving under 
gravity, i.e. S n = constant is preserved by a local biasing, i.e. S^ remains a constant 
(with a different value) at all orders n. On the contrary, as was shown in p]] |, a 
cooperative galaxy formation scheme can be severely constrained by the skewness it 
would predict. 
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Finally, let us mention that a similar scaling of the moments with the second one 
has been found in the APM survey ^5|, which is indeed expected || theoretically. 

5. ZELDOVICH APPROXIMATION 

As was shown above, a lot can be accomplished by using the standard Eulerian 
perturbative approach to gravitational instability. However, perturbative equations 
are often easier to integrate when expressed in Lagrangian coordinates. In our 

below). Another 



7.2. Moreover, 



experience, this happened at second order for fi ^ 1 (see §6.1 
example is the redshift space distortion for (S 3 ), discussed in 
(J 2 ) grows with time, at any fixed order in perturbation theory, the Lagrangian 
approach is likely to remain valid longer than the Eulerian approach. This is so because 
the requirement of small Lagrangian displacements and gradients is weaker than the 
requirement (<5 2 ) -C 1. This idea, which motivated Zel'dovich |Q approximation 
(which is simply the extrapolation of the first order Lagrangian solution in a regime 
where its validity is not mathematically warranted), remains valid at higher orders. 

In the following, I first introduce Zel'dovich approximation by using previous 
results concerning the linear Eulerian theory. I then discuss an application to the 
growth of galactic spin (tidal torque theory), as well as the idea of "Pancakes". 

In a second step, I outline the general derivation of the perturbative (Lagrangian) 
solutions, and focus on the redshift distortion problem as an application. I also give 
examples of comparison between Eulerian and Lagrangian theory when (<5 2 ) ~ 1, i.e 
when Zeldovich and higher order approximations are used as approximations to the 
real dynamics, in a regime where there validity is hard to assess analytically. 

5.1. Derivation from Eulerian Theory 

We found above that in linear Eulerian theory the evolution of the growing mode is 
self-similar, with 

5(x, t) = D(r)e(x), </>(x, r) = — &(x); V 2 ^ = 47rGp b a 3 e(x), 

a 

the subscript i referring to an "initial" state. Note that in an Einstein-de Sitter 
Universe, D{t) cx a, and thus the potential is (linearly) conserved. In addition, the 
linearized Euler equation using the conformal time r, ov + av = —aVcj) = —DV/fii 
(with v = x) yields v = -(l/a) f DdrVfa = -(1/-D) J DdrVcf) (i.e. v is parallel to 
V0). 

We can thus use these expressions to get the linear trajectory of a particle as 
a displacement field, ?t?(t,q) = x q, applied to its initial position, or Lagrangian 
coordinate, q, 



x = q+ M>(t, q) with * 



dr . 
— / DoIt 



a 



(25) 



By using the continuity equation written as aD = aD = AnGpi,a D, we obtain 

*(r,q) = 6(r)* i , with b(r) = - 3 ,^ = V0 t , (26) 

the velocity field being given by v = —ab^i. At the linear order then, particles just 
go straight (in comoving coordinates) in the direction set by their initial velocity. 
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Zeldovich' approximation is to use these linear solutions to extrapolate fluid 
element trajectories well beyond the validity range of the linear Eulerian theory 
(5«1). 

5.2. Galactic Spin Origin 

Zeldovich approximation allows to see quite easily how the tidal field of an 
inhomogeneous distribution couples to the quadrupole of the mass distribution of 
a proto-object to induce a growth of its angular momentum M. At early times, 



M = / d 3 q pa 3 (ax. - ax) A v/a, 
Jv 



where the integral is taken over the initial (lagrangian) volume of the perturbation V, 
of center of mass x = J y d 3 qx./V . At the lowest order, p = p b and v = —abVfa, and 
thus 

M = -p b a 4 b / d 3 q (q - q) A V</> 4 = -p b a 4 b / c^(q - q) A dS. 
Jv Js 

The second expression follows by integration par part (Gauss' divergence theorem), 
and the integral is thus taken over the boundary surface of the lagrangian volume of 
the future object. If this volume is spherical (or its surface S is an equipotential) , 
then M remains zero. Provided the initial potential is smooth enough, we can take 
advantage of a Taylor expansion of the acceleration around the center of mass, 

| q = |q + (q - q) • VV^i |q + O [(q - q) 2 ] , 

to relate the angular momentum to the inertia tensor of the proto-object, I 1 and the 
(initial) tide tensor at the center of mass position, T defined in (@) , 



Mi = -abetjkTjilik, with X ik = (qi - qi)(qk - qk)pba 3 d 3 q 



v 

e standing for the totally antisymmetric tensor of order 3. 

In the Einstein-de Sitter case, the linear growth rate of the angular momentum is 
thus 

M cx t 3 cx t. 

Of course, this description holds only during the early collapse phase (~ while the 
proto-object is still expanding). Later one, particle turn over, and there is no chance 
that their trajectories may be described by a simple ballistic approximation. Latter 
phases have been mostly studied by numerical simulations. Still, the simple exercise 
above shows how Zeldovich approximation clearly unveils the basic principle of the 
tidal torque theory of the galactic spin origin, i.e. that the tidal field of neighboring 
density perturbation acts on the irregular shape of the proto-object. 

5.3. Density Contrast & Pancakes 

The density contrast simply obtains from the trajectories by requiring that mass be 
conserved, i.e. 

p(r)d 3 x = p(Ti)d 3 q 1 + 6 = -, 
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J standing for the Jacobian of the transformation from q to x: 
9x 



J 



dq 



\detT>\ 



The deformation tensor, T> may be written as the identity matrix plus the tidal tensor 
T = \dty/dq\, which is directly given by the curvature of the initial potential field (in 
Zeldovich approximation, T = bW^i/ (^nGpba?)). 

Let us call Ai > A2 > A3 the local eigenvalues of the (initial) tidal field. The 
density contrast may then be written as 



[1 - 6(t)AJ [1 - i(r)A 2 ] [1 - 6(t)A,] 



which shows that, if this approximation remains approximately valid till the first 
crossing of trajectories, one expect a singularity, a (locally) planar collapse to infinite 
density along the axis of the largest positive (if any) eigenvalue, Ai, when 6(r)Ai — > 1. 
If no eigenvalue is positive, this is a developing underdense region, while in the rarer 
case when two eigenvalues are approximately equal, Ai = A2 ^ A3, there is a collapse 
to a filament. The case Ai = A2 = A3 leads to a spherical collapse. 

For a gaussian field, Doroshkevich found the probability distribution of the 
eigenvalues, and in particular that P(Xi > 0, A2 < 0, A3 < 0) = 0.42 ( = P(Xi < 
0, A2 > 0, A3 > ). It shows that indeed one generically expects "pancakes" around 
local maxima of Ai. This is to be contrasted with the linear Eulerian theory whose 
extrapolation would lead to predicting a collapse around the initial maxima of the 
density field (since S oc Si), i.e. around local maxima of Ai + A2 + A3 whose collapse 
is rather spherical. 

The "pancake" theory arose in the context of a "neutrino" dominated dark matter, 
where there is a natural small scale cut-off in the power spectrum of perturbations due 
to their early phase of free-streaming. Then the potential field is quite smooth and it is 
not too surprising that Zeldovich approximation provides a quite accurate description 
of the large scale structures which appear in numerical simulations. Indeed, we only 
have to describe an essentially laminar flow. On the other hand, in a hierarchical 
scenario where small scale structures collapse first, these small scale potential wells 
will undoubtedly curve the particles trajectories and Zeldovich approximation would 
not describe at all the individual particle trajectories. We shall see later that this 
approximation might still be useful, but in an average sense. 



6. LAGRANGIAN PERTURBATION THEORY 

In the previous section, we made contact with a lagrangian description of the dynamics 
by introducing Zeldovich approximation as a consequence of the Eulerian perturbation 
theory at the linear order. In the following, we rather describe the steps involved in 
a systematic pcrturbative approach of Lagrangian theory, as initiated by Moutarde 
et al 1991, |3?| (see also the contribution of T. Buchert in this volume). We then 
proceed with an application to these results to the computing of the skewness in 
redshift space. We shall later come back to using those rigorous perturbative solution 
as approximations to the full non-linear dynamics, thereby extending Zeldovich ansatz. 
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6.1. Perturbative equations 

The Lagrangian perturbative approach introduced by Moutarde et al (1991) proceeds 
in quite a parallel fashion to the Eulcrian one reviewed earlier. As before, we use 
the Newtonian approximation with comoving coordinates (see the contribution of 
S. Matarese in this volume for the full general relativistic approach). Following 
Doroshkevich et al 1973 |l8|| , we replace the standard cosmological time by a new 
time r defined by 

dr oc a~ 2 dt. (27) 

Beware that we kept the same symbol here than above for the conformal time (which 
is defined instead by adr oc dt\). In that case the motion and field equations now read 

x = -V,*, A x $ = p( r )8, (28) 

with no x term. In this equation, dots denote derivatives with respect to t, 
and /3(r) = 47rGa(pa 3 ). By choosing properly the proportionality constant in the 
definition ( p7| ) when, respectively, f2 is equal, smaller or greater than 1, j3 is simply 
given by 

f T 2 + fc(r>)' 

with k(n = 1) = 0, k(n < 1) = -1, and fc(fi > 1) = +1. While r = i" 1 / 3 , when 
Q = 1, one has r = |1 — $^2 1 1 / 2 otherwise. 

Now we wish to follow the particle trajectories instead of the density contrast, by 
using the mapping 

x(r) =q+*(r,q) (29) 

The Jacobian of the transformation from x to q, permits as before to express the 
requirement of mass conservation simply as 

p{x)Jd 3 q = p(q)d 3 q, 

i.e., S = J -1 — 1 (note that p(x) is in fact a shorthand for p(q, t)). By taking the 
divergence of the equation of motion ( |2§| ) , one obtains the equivalent of the Eulerian 
equation (Eh: 

J(r,q)V x x = /3(r)[J(T,q)-l]. (30) 

Of course the addition of any divergence-free displacement field to a solution of 
the previous equation will also be a solution. In the following, we remove this 
indeterminacy by restricting our attention to potential movements, which must satisfy 

V K xx = 0. (31) 

The main reason to restrict to that case is that vortical perturbations linearly decay, 
a consequence of the conservation of angular momentum in an expanding universe. 
Thus one might consider that the solutions will apply anyway, even if vorticity is 
initially present, because at later times it will decay away. In the same spirit, we 
shall mainly focus on growing mode solution (see |l^, |l3) for the cases of rotational 
perturbations in Lagrangian space but not in Eulerian space, and also for the effect 
of decaying modes). 
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The final equation to solve obtains by rewriting the divergence of the acceleration 
r = x explicitly as a function of q 



(32) 



where the Ay are the cofactors of the Jacobian, and the partial derivatives denoted 
by Latin letter are taken with respect to the q coordinate (e.g., , Tij — ^f). 

As in the Eulerian case, perturbative solutions are obtained by means of an 
iterative procedure. But this time, the expansion concerns the particles displacement 
field itself, and we write it as 



* = e* (1 > + e 2 ^ + e 3 *^ + 0(e 4 ). 
The determinant of the Jacobian is then similarly expanded as 
j =1 + £ jW +£ 2j(i) +£ 3jW +0 (s i ) 



(33) 



(34) 



+ £ 3 (K& + £(3) + M (3)) + ( £ 4) 

where K^ m \ L (m \ and denote the m-th order part of the (invariant) scalars 



K 
L 



M = D = det[*ij- 



In other words, 



£(2) 
£(3) 



= v • * (m) = V * (m 

>(2) lT/ (l) _ 



) 



^(2)^(1) 



(35) 



det[*| 1}l 



The equation to solve then follows by replacing the expansions (|33|-p5[) of the 
displacement field and of the Jacobian in equations ( p0| ) and (|32|). We now proceed 
to an iterative solution of this equation. 



6.2. Low order Solutions 

I give below only the solutions for the first two orders, when A — and CI is either unity 
or smaller than one. This is all we need for the application described afterwards. A 
more complete set in solution can be found in Bouchet et al 1995 [|l0| who give for the 
first three orders the solutions for arbitrary values of VL when there is no cosmological 
constant, as well as numerical solutions for flat universes when A ^ 0. 

6.2.1. First Order: keeping only the terms of order e, one gets 

K (1) - (3K (1) = 0. (36) 

Thus i^W(r,q) can be factorized into a temporal and spatial part, K^'(t, q) = 
gi(r)K^ (n, q), where if^^r^q) is the divergence of the initial displacement field, 
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9i( T i) i s assumed to be unity, and g\ = k a gi a + k b gi b . For = 1, the linear growth 
rate g\ is given by g\ a = t~ 2 , and gu = r 3 . When < 1, 

.9 la = l + 3(T 2 -l)(l + rS), fflfc = r(r 2 -l), S = \ In ^ . 

The closed case is given in Q. 

For a potential movement, one thus recovers Zeldovich solution, ^^'(t, q) = 
i?i(t)\I'( 1 ) (q), as in (p5|), if we define (for all m) 

* (m) (q) = * (m) (^,q)- 

Also, note that (5 (1) (r) cx gi(r) (with initially <5 4 (1) = -V • * (1) (q)), i.e., , the Eulerian 
linear behavior is recovered, since g\ = D. 

The logarithmic derivative of the growth factor, fi = ^j-^-, is useful to describe 
comoving peculiar velocities and therefore redshift distortions (see below). Near 0=1, 
a limited expansion of the solution shows that /i x O 4 / 7 where x means "behaves 
asymptotically as". A better analytical fit for /i in the range 0.1 < O < 1 is given by 
fx « O 3 / 5 , as was originally proposed in fl42| . 

6.2.2. Second Order: collecting terms of 0{e 2 ) in equations ( |30| ) and (^), we 

find 

K^-f3K^ = -l3glL^{n^). (37) 

Thus, the solution is again separable, i^( 2 ^(r,q) = g 2 {r)K^ (r^, q). The spatial part 
is given by K^(Ti,c[) — L^(Ti,c[), and the growth factor is 

.92 = k 2 a g 2a + 2k a k b g 3b + fc 2 g 2c + Lgia + hgib ■ (38) 

For O = 1, we simply find g 2a = -%gf a , gib = |.9iaffib, g 2c = ~\gu, while for < 1, 

.92a = l-^(r 2 -l) [r+(r 2 -l)S] 2 , 

g 2b = J(t 2 -1)[t 3 + (t 2 -1)S], (39) 
92c= -\ir 2 -lf- 

If only a growing mode is initially present (k b = = /&), which we assume, one 
must have g 2 /gi — > when r — > oo, i.e., , when O — > 1. A limited expansion of 
equation ( p9[) then shows that ^ a must be equal to — |, and this yields the physically 
relevant solution 

.9 2 = -^-^r 2 -l){l+r5+i[r+(r 2 -l)5] 2 }, (40) 

which behaves near O = 1 as 

g^- 3 -n-^gl (41) 

while g 2 jg\ — > —1/2 when O approaches 0. The O > 1 solutions are in IP]]! . The 
logarithmic derivative of the second order growth rate, f 2 = behaves as 

near = 1. Alternatively, it can be somewhat better approximated 
by 2 4 / 7 for 0.1 <, O <, 1. 

The third order solution is also know and was obtained in jlO| . It is also separable 
in the = 1 case. But when O ^ 1, it must be written as a sum of 2 separable pieces. 
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7. APPLICATIONS 



Here I show how moments of the PDF can be computed in Lagrangian perturbation 
theory, and recover some of the result reviewed above. In particular I give the skewness 
for rather general non-gaussian initial conditions. Then I show how the redshift space 
distortion effect on skewness can be computed. Finally, I discuss using the exact 
perturbative solutions as approximations to the real dynamics. 



7.1. Moments in the Lagrangian approach 

Since the unsmoothed density contrast is given by 8 = J -1 — 1, we have 

<<Hx)) x = ((J- 1 - i)") x = ((J- 1 - i)" J) q ■ 

The variance and skewness are then given up to the fourth order by 

(S 2 ) = e 2 (JW 2 ) + e 3 (2JM j( 2 ) - jW 3 ) + 

e 4 (JM 4 - 3J« 2 J< 2 > + J( 2 ) 2 + 2JM J0)> + 0(e 5 ) 
(6 3 ) = -e 3 (J« 3 ) + e 4 (2J« 4 - 3J« 2 J( 2 )) + 0(e 5 ), 

where all averages on the displacement field are taken with respect to the Lagrangian 
unperturbed coordinate q. 

Now consider the case of the initial density fields which can be generated by a 
displacement field 'S'W which three components are independent, with the same (but 
arbitrary) statistical law to insure homogeneity and isotropy. While not fully general, 
the generated density field Si = eS^ 1 ' encompasses a rather wide class of non-gaussian 
Initial conditions. We note a\ the variance of any component i of the (linear) gradient 
field 



and S its third moment S = /\Jj^ 3 \, and K its reduced fourth moment K = 

(*g 4 )"3af. We thus have < J« 3 ) = 3S, < J« 4 ) = 3K + 27a 4 = 3K + 3 ( J« 2 f . 

The other term in (<5 3 ) involves the product j' 1 ) 2 which can readily be estimated 
since [after (^5|) and (]35|)] we have 

It follows by development that ( J (1) 2 J (2) ) = 6(1 + g 2 jg 2 )a\. We thus have 

-35 + e gi 6(K + 3(2-g 2 /g 2 )at) 
e gi at(3-3sS/o- 2 : 

„^+ 4 - 2. 92 /. 9 2 + 2 ^-^ 

For an initially Gaussian field with 5 = if = 0, we get the simple result 

S 3 = 4-2g 2 /g 2 + 0(e 2 ), (44) 

whose first term corresponds to the pure Zeldovich approximation and had been found 
by Grinstein and Wise |27t| . The value of the ratio of growing modes may be obtained 



S "- ~ T 4 (3-3 £ 5/a 2 ) 2 ■ N "' 

J±-+A-2g 2 /g 2 + 2 ^l~ u +0(e). (43) 
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from the exact solutions given before, but it is hard to imagine practical cases when 
our approximation ( |4l] ) might not be sufficient. We thus have the handy and quite 
accurate formula for a Gaussian initial field with 0.1 Ss fl ^ 2 (Bouchet et al 1992) 

28 + 6 Q- 2 / 63 
- 7 , 

which generalizes the S3 = 34/7 found by Peebles (1976) in the £1 = 1 case. The more 



general formula (43) gives the time-evolution of the skewness factor S3 in spatially 
flat models for the class of non-Gaussian initial conditions which may be generated by 
independent displacement fields along three axes. The general formula for the Q = 1 
case can be found in . 



7.2. Redshift Space Skewness 

In redshift space, the appearance of structures is distorted by peculiar velocities. At 
"small" scales, this leads to the "finger of god" effect: the clusters are elongated along 
the line-of-sight due to their internal velocity dispersion. This is an intrinsically non- 
linear effect, and we shall not be concerned with it. At "large" scales, the effect is 
reversed: the coherent inflow leads to a density contrast increase parallel to the line-of 
sight. Indeed, foreground galaxies appear further than they are, while those in the 
back look closer, both being apparently closer to the accreting structure (Sargent & 
Turner 1977; LSS, §76; Kaiser 1987). 

Kaiser (1987) estimated this redshift effect, in the large sample limit, on 
the direction averaged correlation function, and found (<5 2 ) = (1 + 2/3 O 06 + 
1/5 fi 1 ' 2 ) (<5 2 ), where the superscript z correspond to a redshift space measurements. 
For fl = 1, {5 2 ) z = 28/15 (5 2 ). We now turn to a Lagrangian analysis of this effect. 




-3 -2-10 1 -3 -2 -1 

n n 



Figure 5. See legend in text. Courtesy Hivon et al 1993) 



Let us consider the case of spherical coordinates, when distances to the observer 
would be estimated by means of redshift measurements. And let us now denote redshift 
space measurements by the superscript z. The redshift space comoving position x z of 
a particle located in r(q) = ax(q) is x z = r/(aH) (with H = a/a). The real space 
perturbative expansion (|2^) is then replaced by 



x z =q + [l + /i(t)] 5 i(i)* (1) (q) 

+ [l + / 2 (i)] 5 2W* (2) (q)+0(e 3 ), 



(45) 
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where we have explicitly used the separability of = g 1 (t) 1 4'( 1 )(q) and \&( 2 ) = 

c/2 (t)"*^ 2 ^ (q) . In the limit of an infinitely remote observer, say along the T3-axis, 
we can approximate spherical coordinates by Cartesian ones. The observed density 
contrast S z in comoving coordinates is then simply S z (xi, Xi, x$). All we have to do, 
thus, is to replace everywhere in the calculation of S3 the displacement field along the 
third axis 

*< m > =<,„,*<"*> by (l + f m )g m ^ 

for m = 1 and 2. We have for instance, ( jW 2 ) = [2 + (1 + A) 2 ] a\. 

This shows that in our so-called "infinite observer limit", we have {S 2 ) z — 
(1 + 2/3 /1 + l/3/ 2 )((5 2 ) which slightly differs from Kaiser's calculation (who did 
a calculation in spherical coordinates instead of our rectangular ones). Indeed, we 
find that the redshift space variance is boosted by a factor 30/15 for fi = 1, instead 
of his value of 28/15. 

The calculation of the terms involved of the skewness calculation is fairly 
straightforward, with the exception of a piece £ which is small and can be bounded. 
The result for our fairly general non-gaussian model is 

5 *_ 6 (2+(l + /i) 3 )6' 
3 e(2+(l + / 1 ) 2 ) 2 af 



2 (2+(l + / 1 ) 3 ) 2 S 
(2 + (l + / 1 ) 2 ) 3 a6 
2 (2 + (l + / 1 ) 4 )i^ 
(2+a + A) 2 ) 2 af 
6 1 + 2(1 + A) 2 + (92/ 'gj) [3 + 2fr +f 2 + hh£\ 
[2 + (l + /i) 2 ] 2 

with 1 > E > 0. Of course, we recover the real space result (Q) if we set /1 = / 2 = 0. 
On the other hand, if fi = 1, we have /1 = 1 = A/2 (and gilg\ = — 3/7), which yields 
for Gaussian initial conditions S 3 = (35 + £)/7 while, for n = 0.1, S 3 « (34.5 + 0.4f )/7 
(to be compared with the value of 34/7 in real space). 

The formula ( ^6|) obtained above applies only in the limit of large volumes (like 
Kaiser's result), since we have taken the limit of an infinitely remote observer. A full 
calculation along the previous lines, but in spherical coordinates, and including the 
effect of smoothing, is done in p8|| . In any case, equation ( |46| ) clearly shows that the 
ratio 5*3 is, for Gaussian initial conditions, nearly independent of the value of and 
is barely affected by redshift space distortions. It is interesting to note, though, that 
in the non-Gaussian case, the distortion might be rather large, for large enough S or 
K. 
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7.3. Approach of the Non- Linear Regime 

So far, we have mainly considered rigorous uses of the Lagrangian perturbative 
approach, for instance the derivation of a second-order quantity, the skewness of a 
PDF, with the help of second order perturbation theory. 

Now, we examine to what extent the second-order theory brings improvement to 
Zeldovich approximation, when both are used as approximations to the real dynamics, 
i.e. outside of their rigorous validity range. 
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The first and second order solutions in Eulerian and Lagrangian perturbation 
theory were compared in |l0[ to spherically symmetric cases whose evolution is 
analytically known. For instance, the results of those approximations were checked for 
the density and the divergence of the velocity field in the spherical top-hat case, when 
its amplitude is varied. Other comparisons were made also for given smooth profiles. 

Of course, direct comparisons with numerical simulations (e.g. Ref. pp| or Buchert, 
this volume, and references therein) were also performed. It turns out that the 
Lagrangian approximations perform considerably better on a given non-linear scale I 
if they are applied to initial conditions smoothed on that scale. This avoids amplifying 
the small scale "noise" whose structure cannot be followed anyway. 

These comparisons lead to the following approximate ranking (at least for 
moderate final density contrasts ~ 1): 



Here the signs ">" and "~" mean respectively "more accurate than" and "of 
comparable accuracy to" . For relatively large final density contrasts, the Eulerian 
approach becomes particularly inefficient, except for the velocity field, for which 
it tends however to be less accurate than the Lagrangian one. The second order 
Lagrangian approach gives, for moderate final S, an excellent approximation of the 
density contrast and the velocity field. Its seems to be able to reproduce density 
contrasts as large as ten. 

There is hope that these lagrangian approach to the non-linear regime might lead 
us to a better understanding of this most poorly known of all regimes, the transition 
between the weakly and strongly non-linear realm. It has been recently an area of 
considerable activity and we shall have to wait till the dust settles. 

8. CONCLUSIONS 

Perturbation theory already has a long and distinguished history in Cosmology, but 
for many years the attention focused mainly on the linear terms, e.g., for early universe 
calculations of the power spectra of light and matter at recombination. But for the 
latter evolution of large scale structures, perturbation theory went out of fashion in the 
era which saw the advent of massive computer simulations of their formation. These 
allowed investigating strongly non- linear scales, the only ones that could be carefully 
studied statistically in the galaxy surveys of the time. Another reason of this purgatory 
period of the theory may have been the perception that any early signature of the 
weakly non-linear phase would have been totally erased by the following strongly non- 
linear phase, i.e., that it was irrelevant. In any case, the situation has now dramatically 
changed, galaxy catalogs encompass ever increasing scales where density contrasts are 
weak; they display remarkable moments hierarchy on large scales, with a striking 
similarity to that observed on smaller scales. Meanwhile realistic calculations have 
become recently available (e.g., including the smoothing inevitable in any realistic 
measurement, redshift distortions, etc), and last but not least, comparisons with 



density contrast : 

Lagrangian second order 
Eulerian second order 



> Zeldovich > 

> Eulerian linear theory, 



Velocity field : 

Lagrangian second order 
Zeldovich 



Eulerian second order > 
Eulerian linear theory. 
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numerical simulations showed the predictive power of the theory. In a few years the 
beauty awakened, and this area of research is now striving. These notes just surveyed 
some of the simpler analyses, and new results are posted nearly every week. Of course, 
much remain to be done, and maybe the most important question to be answered is 
why does these perturbative approaches work so well? 
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